Electro-osmotic flows under nanoconfinement: a self-consistent approach 
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We introduce a theoretical and numerical method to investigate the properties of electro-osmotic 
flows under conditions of extreme confinement. The present approach, aiming to provide a simple 
modeling of electrolyte solutions described as ternary mixtures, which comprises two ionic species 
and a third uncharged component, is an extension of our recent work on binary neutral mixtures. 

t-H The approach, which combines elements of kinetic theory, density functional theory with Lattice- 

t-H Boltzmann algorithms, is microscopic and self-consistent and does not require the use of constitutive 

equations to determine the fluxes. Numerical solutions are obtained by solving the resulting coupled 
equations for the one-particle phase-space distributions of the species by means of a Lattice Boltz- 

,— i mann discretization procedure. Results are given for the microscopic density and velocity profiles 

^ and for the volumetric and charge flow. 
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The issue of transport of ionic solutions confined in nanometer scale channels has been a concern of physicists, 
biologists and chemists for a long time due to their importance in technological, biological and industrial applications. 
With advances in nanofabrication techniques the study of fluid dynamics in atomically sized systems is now accessible 
r>2 and has become an active field of research [TJ[2]. 

Electro-osmotic flow is an important phenomenon occurring in electrolytes when an electric field is applied along 
a channel and drags the electric double layer (EDL) that spontaneously forms next to the charged walls. Due to the 
existence of Debye screening, the resulting velocity profile of the fluid strongly departs from the classical Poiseuille 
velocity distribution, being plug-like, with the consequence of greatly reducing the viscous dissipation near the channel 
walls. While the understanding of such phenomena in microchannels is almost complete, in nanochannels the picture 
is less clear, because the density, charge and velocity profiles can be highly non-uniform at such reduced scales [3 . 

The basis to understand the transport of electrolytes in narrow capillaries was laid down by Smoluchowski about 
a century ago [¥|. This theory is a mesoscopic treatment based on the combination of two continuous approaches: 
the formation of a non-uniform spatial charge near the charged boundaries, the EDL, is described via the Poisson- 
Boltzmann equation whereas the mass flow induced by an applied electric field is determined through the Navier-Stokes 
equation. 

In spite of their crude underlying assumptions, the Smoluchowski approach and its modern evolutions capture most 

r-o of the physics associated with ion transport in microfluidic systems and represent a extremely useful workhorse in 

electrokinetic theory. On the scales characterizing the microchannels (1/im to 100/im), the equations of hydrodynamics 

together with the Poisson-Boltzmann equation describing the electric charge distribution provide a satisfactory picture 

of the interplay between diffusion of the various species, Coulomb interactions and mass transport. The coupled 

£ — equations must be solved with the spatially nonuniform boundary conditions needed to describe the presence of 

charges at the surface of the channel. Nevertheless, the validity of such an approach is questionable when the pore 

size becomes comparable either with the size of the constituting molecules or with the typical length characterizing 

the behavior of electrolytes, the Debye length A D = l/k D . 

^ A fully microscopic investigation of the system via atomistic simulation, such as Molecular Dynamics, is feasible, 

but while it provides an accurate description of the short-time, small-scale behavior, due to the intrinsic noise present 

in simulations, it precludes to access large scales where macroscopic laws emerge and conditions of sub-millimolar 

concentrations have to be considered [5]. 

About forty years ago, Gross and Osterle (GO) formulated the very popular space-charge model where the elec- 
trostatic potential is determined by means of the Poisson-Boltzmann equation, the diffusion of the charges by the 
Nernst-Planck equation and the velocity for the laminar flow by the Navier-Stokes equation |B]. The approach pro- 
posed by Nilson and Griffiths |7| partially goes beyond the ideal gas treatment of the ions of GO and introduces an 
explicit treatment of the solvent, by employing tools of density functional theory to determine the density of the three 
species, but still considers the mean local velocity of the fluid at the macroscopic Navier-Stokes level and the trans- 
port coefficients are not not determined self-consistently and microscopically, but introduced as phenomenological 
parameters. 

In the present letter, we show that it is possible to determine within a unified framework thermodynamic and 
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structural properties and dissipative forces starting from a microscopic description of the interactions among the 
particles in a simple model of ionic solution. For this purpose, a natural candidate to bridge the gap between 
macroscopic and microscopic approaches which also treats on equal footing density and velocity inhomogeneities is 
the modified Enskog-Boltzmann equation [llj . There, the complete information on the statistical description of a 
fluid is assumed to be encapsulated in the one-particle phase space distribution function of each species present in the 
fluid. From the knowledge of the distribution functions one computes the mean values of the physical observables of 
interest, such as charge, mass, velocity and temperature distributions. Within this framework, hard-spheres represent 
a central reference to handle memory- free collisional rules by means of a modified Enskog kernel [12j . 

II. GOVERNING EQUATIONS 

Let us consider for the sake of simplicity a ternary mixture composed by hard spheres with a common diameter 
a and equal masses m. The model can be easily extended to include species of different diameters and masses. We 
associate positive and negative positive charges, z±, to the species labelled + and — , respectively, whereas the third 
component with label carries zero charge. The interactions between point charges and between these and the surface 
charges at the boundaries occur in a space characterized by a uniform dielectric constant e [T5] . 

The electrostatic potential <p(r, t) acting on the charges at point r is given by the contribution determined by the 
ionic charge distribution p(r, t) via the Poisson equation: 

V V (M) = _£M (1) 

and to the "external" charges sitting on the electrodes. No retardation effects are included in our treatment. 

To describe dynamically the system we associate to each species the one-particle distribution function f a (r,v,t) of 
each component a, with a = +, — , 0, whose evolution is governed by the following set of integro-differential equations: 

8 V a (r) 8 P7 8 

|>(r,v,i) + v • Vf (r,v,i) + — ^ ■ £/ a (r,v,i) - -^V^(r) ■ £-/ a (r,v,i) = 
at m av m av 

-c[r(r,v,t)-V^(r,v,i)]+/3* Q (r,t)-(v-u(r,i))^(r,v,i). 

(2) 

The number density n a (r,t) and the average velocity u" of species a are obtained by multiplying by 1 and v, 
respectively, the distribution /" and integrating w.r.t. v. The average packing fraction is £ = 7r/6cr 3 ^ Q n a and the 
local barycentric velocity of the composite fluid is: 

, , T n a (r,t)u a (r,t) 
u(r,t)= ^° V ' H-A (3) 

The temperature ksT = 1//3, with ks being the Boltzmann constant, is assumed to be uniform throughout the 
system and each species has thermal velocity vt = \JksT /m. In addition, F" and <fr Q represent the sum of external 
and internal forces respectively of non-electrostatic nature acting on species a. It should be noted that electrostatic 
interactions are introduced at mean-field Vlasov level and due to the nature of the approximation involved do not 
contribute to the dissipative character of the fluid, while they contribute to the chemical potential |8J. As discussed 
in ref. [2], adding a Coulomb correlation to the local chemical potential does not significantly alter the ionic and 
solvent densities as compared to the mean field approximation for electrostatics. 

The l.h.s. of eq. (pi encodes the effect of streaming on species a while the r.h.s. contains the sum of two terms, a 
first one being a modified Bhatnagar-Gross-Krook relaxation kernel towards a Maxwellian |5J|TD], with 

r(r^t)=n^, t)j ^e^J^^) (4) 

and 
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where the collision frequency v relates to the ideal contribution to viscosity. To summarize the meaning of eqs. (|2|5|, 
the viscosity of the system must contain two contributions: the ideal viscosity, that we represent by the BGK term, 



and the hard-sphere viscosity, that we represent by the last term of eq. (pi. This second contribution is strongly 
modulated by the fluctuations in density and velocity, whereas the BGK term is not. As far as the meaning of the 
baricentric velocity at atomic scale is concerned, one can consult the work of Travis and Gubbins [3|. 

The second term in the r.h.s. of eq. (p| is a collisional contribution. Since all interactions, but the hard-core 
repulsion are treated within the mean-field approximation, the form of <fr a solely depends on the properties of the 
underlying hard-sphere mixture. Therefore, we can use the results already established for neutral systems and derived 
from a simplified treatment of the Revised Enskog Theory (RET) [TT] due to Santos et al. |15j . 

Deriving the mesoscopic equations for the model above is a relatively simple task and can be achieved by using 
the same methods presented in refs. (21 HE] One obtains an equation for the charge distribution, an equation for the 
velocity distribution, assumed to be the same for all components. The continuity equation reads 

^n Q (r,t) + V • (>(r,i)u(r,i))+V ■ (n a (r, t) (u a (r , t) - u(r,t))= 0, (6) 

The propagation of momentum of the single species obeys the following equation 

^K(r,^(r,i)] + ^(n Q (r,0<(r,i)u"(r,t)-n Q (r,0(<(r,i)- U .Kr,i))K(r,i)~ Uj (r,t))) = 

d 7rg(r,t) F/(r) g?(r,t) ez a d 

-- J - + - 1 n a (r,i) + - 1 n a (r,t) n a (r,t) — 0(r,i (7) 
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where 

Kj( r ^) =m dw{v t - Ui)(vj - Uj)f a (r,v,t) (8) 

is the kinetic contribution of component a to the pressure tensor. As a result of a detailed microscopic calculation one 
can show that the average local force acting on a particle of species a due to the influence of the remaining particles 
can be subdivided into three different contributions of non electrostatic nature, associated with the excess chemical 
potential, the drag and viscous forces: 

* Q (r, t) = F a - m/ (r, t) + F a ' dra9 (r, t) + F a ' visc (r, t). (9) 

The explicit form for these quantities depends on the microscopic details and for the charged hard-sphere model 
employed expressions can be found in ref. |5]. 

III. RESULTS 

The numerical solution of eq. (p| is obtained by a generalization of the method presented in ref. [16 to the 
case of a ternary mixture, following the ideas outlined in ref. [17 . Similarly to the Lattice Boltzmann method 
of solving the BGK dynamics [191 EU] , our numerical approach spatially discretizes the distribution function over a 
cartesian mesh. The velocity dependence of the distribution is represented according to a truncated Hermite expansion 
supplemented by a Gauss-Hermite quadrature to evaluate the kinetic moments |18j . The end result is a replacement of 
the distribution function with a discretized version, / Q (r, v, t) — > /5(r, t), where the index p labels a set of discretized 
velocities {c p }. In the following, we use the so-called D3Q19 mesh, including 19 discrete speeds related to exchanges 
of particles between zeroth, first and second mesh neighbors. The Hermite expansion introduces a characteristic speed 
that we identify with the thermal speed vt, being the same for each fluid component. 

The numerical accuracy of the simulations is controlled by the truncation level of the Hermite expansion and the 
associated Gauss-Hermite quadratures in velocity space. The collisional contribution in eq. &h and the Maxwellian 
in the BGK-like term are expressed up to second order Hermite components. The number density and fluid velocity 
are obtained as n a = X^»/p ano - n a u a = ^2„c p fp, respectively. Finally, external forces arc implemented according 
to the method suggested in ref. |^U m order to achieve second-order accuracy in time. 

The method employed to solve the ternary system is further complemented with the solution of an auxiliary 
distribution to minimize the effect of parassitic currents in the system. In fact, given the strong density gradients in 
the system, the LBM method has the attitude to develop spurious currents in proximity of confining walls, where the 
density of species has the strongest variations. This limitation is intrinsic to the spatial discretization of the Lattice 
Boltzmann method. To solve this issue we have developed an ad-hoc strategy to circumvent such problem, that will 
be presented elsewhere [22 . The main idea is to use an implicit trapezoidal rule for the time integration together 
with a version of the two-distribution method of He et al. 123 . 



The systems are simulated in three-dimensions in slit geometries of variable width W along the z direction, with 
—W e ff/2 < z < +W e ff/2 where W e ff = W — u is the available width, and being periodic in the x and y directions. 
To resolve accurately the density profiles, we have chosen a hard sphere diameter to be equal to 8 in lattice spacing 
units. At positions z = ±We//, no-slip boundary conditions are imposed for the fluid velocities via the bounce-back 
method applied on the populations of each species [TS] and at this level, we disregard slippage effects together with 
the distinction between the Stern layer and the shear layer. We consider a 1:1 symmetric electrolyte by enforcing 
electroneutrality on the composite saline solution and charged surface system. The Poisson equation is solved numer- 
ically at every iteration of the LB algorithm by a successive over-relaxation (SOR) method. Second-order accurate, 
von Neumann boundary conditions on the electric field are imposed, that is, h ■ d(f>\ z= ±w = — E/e, where E is the 
superficial charge density on the plates. The Bjerrum length, Ib = e 2 /(AireksT), is taken equal to Ib = 1.5c, a value 
compatible with bulk water conditions at ambient temperature, and data are collected for a uniform electric field 
na 3 eE/mr]f d — 0.005 acting in the direction parallel to the walls. 

We remark that the mean free path of a hard sphere system in bulk is A m / = -t= — } 2 . By confining the system 

in a slit of width W, the Knudsen number expresses the ratio of the two characteristic lengths of the problem, that 
is, Kn = -p££ ~ 0.1 for the systems that we consider in this letter. In other words, collisions are frequent enough to 
restore local equilibrium. 

The density and velocity profiles for two values of the surface charge are shown in Fig. [T] The profiles are obtained for 
packing fraction £ = 0.28 and are compared with the case of negligible packing fraction, whose evolution corresponds 
to electrokinetic motion in the absence of steric forces. The reported normalized densities are n a '* = n a /n,Q where 
Uq is the corresponding bulk value, whereas the normalized velocities are u a -* — u a /uq, with uq — eE / '(Vid^t b) and 
rjid being the ideal contribution to the dynamic viscosity of the mixture. 

The data in Fig. ffl demonstrate the strong oscillations in the density in proximity of the walls and the large 
differences with the case of the ideal electrolyte, with the characteristic cusp of the hard-sphere fluid at contact with 
the wall. The oscillations are further modulated by the development of the Debye layer, whose strength increases with 
the surface charge. Analogously, the velocity profiles develop non-trivial shapes that are far from quasi-parabolic, as 
compared to the ideal solution. It is worth noticing both the peaks in the velocities of the charged species close to the 
wall, and the plug-like shape in the inner part of the slit. The inhomogeneities in the velocity profiles denote major 
departures from the ideal solution case with a reduction of the velocities in the bulk region and an overshooting near 
the walls. As the figure reveals, as a general trend the barycentric velocities have a profile that stays closer to that of 
the neutral species than in the ideal case due to a drag effect. 

The competition of the three species arising from drag and viscous forces produces the net flowrates shown in 
Fig. [2j Data are reported for different pore widths and for different values of the Debye length, ranging from thin 
layer (A D <C W e ff) to overlapping Debye layer conditions (A D > W e ff). The normalized flowrates, computed as 
/ dzu* z (x), are compared with the prediction arising from the Debye-Hiickel (DH) level of theory, being 
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(10) 



The data show large departures from the DH prediction, with a strong suppression of flowrates at increasing slit 
width due to steric effects. The results of Fig. [2] demonstrate the non-trivial effects of nanomctric crowding, where 
the flowrates decrease as the packing fraction increases with a strong dependence on the degree of confinement. 

In nanofluidic applications, it is important to evaluate the role of steric interactions in evaluating the balance of 
bulk dominated versus surface dominated charge current. We then compute the electro-osmotic charge current, 
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iff 

Ieo= I (n + (x)-rr(x)K(x) (11) 



which is of convective origin, and compare it to the conductive component, Iohm = I — IeOi where the total current 
is measured as 

fW e ff 

1=1 dz(n + (x)u+(x)-n"(x)u"(x)) (12) 
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Fig. [3] reports the ratio of currents and shows a large dominance of the Ohmic component stemming from steric 
interactions for all values of the slit geometry and for two values of the surface charge considered, as compared to 
the virtualy zero packing fraction case. These results provide important indications on the effect of ionic crowding in 
nanospaces, as much as providing a guide to experiments and theoretical investigations in nanofluidic set-ups. 



IV. CONCLUSIONS 

In summary, the results presented in this letter provide a working framework for applying theoretical analysis and 
simulation to nanofluids. The method delivers unique informations relative to the transport of ionic solutions confined 
at molecular scale. The theory presented in this letter is self-consistent in such a way as to include thermodynamic, 
drag and viscous collisional forces in a unified way and avoiding the use of fudge parameters. The simulation is stable 
and robust, as it delivers the noise-free distribution of the hydrodynamic moments for a wide range of parameters. 
Work is in progress to analyze systematically the several consequences of collisional forces in electrokinetics. 
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Figure 1: Profiles of solvent (black), counterions (red) and coions (green) for 0.28 (solid lines) and virtually zero (dashed 
lines) packing fractions. The left and right panels correspond to normalized density profiles and normalized velocity profiles, 
respectively. The upper and lower panels correspond to surface charge Ecr 2 /e = 0.0064 and E<r 2 /e = 0.064, respectively. The 
data are for a slit with W = 64 and k D W = 5.33 
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Figure 2: Normalized flowrates obtained for varying slit width 144// at fixed fc D = 1/1. 5a (left panel) and varying k D at fixed 
W e ff/a = 8 (right panel). Data in the left panel are collected for different packing fractions (~ 0.0 for circle, 0.14 for square 
and 0.28 for diamond symbols). The Debye-Hiickel prediction in a equivalent channel width W e ff is reported as a dashed 
line. Data in the right panel are collected for Eer 2 /e = 0.0064 and for different packing fractions (~ 0.0 for circle, and 0.28 for 
diamond symbols) and for Ea 2 /e = 0.064 and for different packing fractions (~ 0.0 for square, and 0.28 for triangle symbols). 
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Figure 3: Log-log plot of the ratio of electrophoretic vs electroosmotic components to the ionic current as a function of the 
slit width. Data are obtained for packing fraction ~ 0.0 (open symbols) and 0.28 (filled symbols) and for surface charge 
£cr 2 /e = 0.0064 (circles) and 0.00064 (squares). Inset: the same ratio reported as a function of the surface charge and for 

W/a = 8. 



